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Fuel  cells  allow  for  increased  efficiency  in  power  production  when  compared  to  the  thermodynamically 
limited  efficiencies  of  heat  engines.  In  the  case  of  solid  oxide  fuel  cells,  they  are  also  usable  with  the  fuel 
infrastructure  currently  in  place  (natural  gas).  Although  potentially  transformative,  solid  oxide  fuel  cells 
are  currently  limited  by  engineering  challenges  related  to  operating  temperature  (>600  °C),  durability, 
and  load  following  ability.  For  example,  the  buildup  of  solid  carbon  in  the  stack,  or  coking,  potentially 
limits  one  of  the  most  desirable  aspects  of  solid  oxide  fuel  cells,  which  is  their  robustness  to  fuel  type.  As 
the  working  temperatures  for  SOFCs  continue  to  decrease,  in  order  to  maintain  fuel  robustness,  the  need 
for  control  of  the  inlet  fuel  composition  increases.  This  work  demonstrates  the  use  of  a  model  predictive 
control  algorithm  on  a  solid  oxide  fuel  cell  system  including  reformer,  blowers,  heat  exchanger,  tail  gas 
burner  and  stack.  The  controller  allows  for  load  following  demand  changes  from  the  fuel  cell  and  meets 
those  demand  changes,  while  ensuring  that  the  reformate  composition  is  not  prone  to  solid  carbon 
formation.  The  controller  meets  current  demand  changes  to  within  0.1  A  s_1  while  maintaining 
compositional  limits  on  the  reformate  flow  and  temperature  limits  on  the  stack  and  reformer. 

©  2012  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

The  use  of  solid  oxide  fuel  cell  (SOFC)  systems  in  both  portable 
and  fixed  installation  continues  to  be  the  subject  of  significant 
research.  The  high  temperatures  required  for  electrolyte  conduc¬ 
tivity  remain  an  engineering  challenge.  Current  systems  employing 
yttria  stabilized  zirconia  electrolytes  tend  to  operate  at  approxi¬ 
mately  700—800  °C.  At  these  high  temperatures  internal  reforming 
of  hydrocarbon  based  fuels  is  possible  [1],  However,  solid  carbon 
formation,  or  coking,  is  still  a  threat  to  system  durability  and 
robustness.  Much  current  work  is  focused  on  lowering  the 
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operating  temperatures  of  SOFC  systems.  Lower  temperatures 
result  in  a  more  efficient  electrochemical  process  and  significantly 
lessen  the  costs  of  external  components  such  as  blowers,  heat  sinks, 
and  transport  systems.  Currently  proposed  SOFC  electrolyte  types 
allow  for  operation  at  temperatures  as  low  as  approximately  400  °C 
[2],  At  these  temperatures  external  fuel  reforming  becomes  even 
more  important  to  avoid  coking  in  the  stack,  as  the  level  of  internal 
reforming  is  reduced  [2],  The  need  for  external  reforming  is  even 
greater  when  biogas  is  the  fuel  source  [3],  The  goal  of  this  research 
is  to  examine  if  a  control  algorithm  can  be  designed  to  control  the 
flow  rates  of  fuel  and  air  to  the  stack  such  that  the  demand  current 
is  produced,  coking  does  not  occur  and  the  operational  limits  of  the 
fuel  reformer  are  maintained.  The  system  to  be  controlled  is 
a  potentially  mobile  1.5  kW  SOFC  system.  Catalytic  partial  oxidation 
(CPOX)  reforming  has  been  chosen  as  it  allows  faster  response 
times,  and  the  addition  of  steam  is  not  required.  The  use  of  steam  or 
auto-thermal  reforming,  although  more  efficient,  decreases 
response  time  and  can  be  problematic  for  mobile  systems.  Air  is 
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assumed  to  be  available  and  both  fuel  and  air  are  moved  through 
the  system  using  blowers. 

Control  of  SOFC  systems  is  challenging  due  to  non-linear  system 
dynamics  and  the  existence  of  multiple  temperature  and  compo¬ 
sitional  operating  constraints.  System  components  have  tempera¬ 
ture  based  operating  limits,  and  the  stack  has  fuel  inlet  composition 
and  utilization  constraints,  including  avoiding  solid  carbon 
formation.  Model  predictive  control  (MPC)  is  a  proven  control 
method  for  meeting  constraints,  but  requires  a  model  that  can  be 
optimized  in  real  time.  Nonlinear  models  that  capture  these 
constraints  are  too  complex  to  optimize  in  real  time  within  a  model 
predictive  controller.  The  approach  taken  in  this  work  is  to  create 
a  first  principles  model  that  captures  the  operating  constraints, 
perform  a  system  reduction  to  obtain  a  set  of  linear  parameter 
varying  (LPV)  models,  and  implement  model  predictive  control  on 
the  resulting  LPV  model. 

The  constraints  captured  include  temperatures  of  inlet  and 
outlet  flows  from  all  components,  spatially  varying  temperature  of 
the  cell  tubes  within  the  stack  and  heat  exchanger  and  a  measure  of 
whether  the  fuel  supplied  from  the  reformer  to  the  stack  is  prone  to 
solid  carbon  formation.  These  constraints  provide  the  motivation 
for  a  set  of  identification  operations  carried  out  using  the  first 
principles  model.  This  simulation  produces  input-output  data 
which  is  then  used  to  identify  the  LPV  model.  In  effect,  the  high 
order  non-linear  model  is  used  as  a  simulated  experiment  to 
identify  the  LPV  model.  The  LPV  model  is  actually  a  set  of  linear 
models  which  are  combined  or  blended  using  a  function  of 
a  measured  signal.  A  linear  model  is  calculated  at  each  time  step  for 
which  the  scheduling  variable  is  measured.  This  linear  model  is 
then  used  to  solve  the  MPC  optimization  problem  to  find  the  new 
inputs  for  the  system  when  the  controller  is  implemented. 

The  proposed  method  could  easily  be  implemented  for  any  fuel 
reformer  type  or  stack  geometry  for  which  an  accurate  first  prin¬ 
ciples  model  exists.  Hybrid  identification  schemes  are  also  possible, 
where  a  specific  fuel  reformer  geometry  or  blower  is  modeled  and 
used  to  produce  input/output  data.  The  reformer  model  outputs 
can  then  be  experimentally  simulated  and  applied  to  a  physical 
system.  The  resulting  input/output  combination  can  then  be  used 
to  identify  a  system  wide  model  for  use  in  model  based  control. 
Thermal  effects  between  modeled  and  experimental  components 
will  not  be  captured  for  such  a  hybrid  identification  scheme. 

2.  Fuel  cell  system 

The  first  challenge  in  developing  a  model  to  be  used  for  model 
based  control  is  determining  the  level  of  fidelity.  Extremely  high 
order  computational  fluid  dynamics  (CFD)  models  can  be  created  to 
match  very  specific  geometries  and  such  models  capture  both  the 
spatial  and  time  varying  aspects  of  the  system.  Such  models  are 
tailored  very  specifically  to  the  geometries  involved  and  impose 
a  prohibitive  computational  burden.  As  such  they  are  most  suited  as 
analysis  and  design  tools  rather  than  for  control.  SOFC  systems  can 
also  be  modeled  using  lumped  thermodynamic  models,  enforcing 
mass  and  energy  conservation,  but  ignoring  spatial  dynamics  and 
chemical  kinetics.  Such  models  also  ignore  the  kinematics  at  the 
reformer  catalyst  and  the  cells.  These  effects  are  important  to  the 
operation  of  the  SOFC  system  and  cannot  be  ignored.  Hot  spots  can 
occur  spatially  along  cell  tubes  [4],  and  the  reformate  composition 
is  highly  dependent  on  the  reaction  kinetics  at  the  catalyst  [5],  Both 
these  factors  impose  operating  constraints  on  the  system.  The 
model  utilized  within  the  controller  has  to  suitably  capture  such 
dynamics  while  still  being  fast  enough  that  solving  an  optimal 
control  problem  in  real  time  is  feasible. 

In  order  to  capture  the  dynamics  of  the  system  as  closely  as 
possible,  we  initially  develop  a  non-linear  first  principles  model  of 


the  system.  This  non-linear  system  model  consists  of  the  following 
(all  components  are  sized  for  a  kilowatt  order  system): 

•  Blowers:  one  to  provide  air  to  the  fuel  reformer,  and  one  to 
provide  air  to  the  stack. 

•  Fuel  reformer:  a  catalytic  partial  oxidation  reactor  which 
produces  the  hydrogen  rich  gas  used  as  a  fuel  for  the  stack. 

•  Fuel  cell  stack:  a  collection  of  tubular  solid  oxide  fuel  cells 
which  produces  the  desired  current. 

•  Tail  gas  burner:  simple  burner  to  combust  any  remaining  fuel 
products  in  the  stack  exhaust  in  order  to  help  pre-heat  air. 

•  Heat  exchanger:  a  counter-flow  tubular  heat  exchanger  to  capture 
heat  from  the  stack  exhaust  and  preheat  the  stack  inlet  air. 

The  blowers  are  modeled  via  first  principles  as  in  [6],  The  blower 
model  is  based  off  a  commercially  available  blower,  an  EBM 
D1G133-DC13-52,  which  has  been  rescaled  using  the  fan  laws  [7], 
The  result  allows  for  the  calculation  of  the  mass  flow  as  a  function 
of  the  motor  speed  and  the  required  pressure  differential.  Manu¬ 
facturer  data  is  used  to  fit  a  function  for  the  motor  speed  and  the 
power  provided  to  the  blower.  No  blower  is  used  for  the  fuel  flow  to 
the  reformer  as  it  assumed  to  be  under  sufficient  pressure. 

The  fuel  reformer  is  considered  as  a  continually  stirred  tank 
reactor  with  surface  chemistry.  A  CSTR  was  used  as  it  is  the  lowest 
dimensional  model  that  can  capture  the  temperature  and  compo¬ 
sition  of  the  reformate.  Although  a  plug  flow  model  would  be  more 
accurate,  only  the  final  exit  composition  and  temperature  are  used. 
The  reformer  has  a  tubular  geometry  with  a  rhodium  catalyst  on 
a  foam  monolith.  The  model  makes  use  of  the  reaction  mechanism 
for  natural  gas  over  rhodium  developed  in  [5].  The  model  itself  has 
been  validated  against  an  experimental  CPOX  reactor  at  the  Colo¬ 
rado  Fuel  Cell  Center  [8],  Reformate  composition  and  temperature 
are  considered  as  functions  of  the  input  fuel  compositions, 
temperatures  and  the  mass  flow  rate  and  temperature  of  the  air 
provided  to  the  reformer.  For  the  purposes  of  identification  and 
control  the  input  fuel  is  biogas  with  a  composition  of  70%  CH4  and 
30%  CO2.  Biogas  was  used  as  an  input  fuel  to  test  the  system  with 
a  high  carbon  content  fuel,  challenging  from  a  coking  standpoint. 
The  biogas  content  was  chosen  as  70%  CH4  and  30%  CO2  since  this 
composition  is  within  the  range  of  that  produced  by  sewage  plants 
and  agriculture.  The  model,  however,  fully  supports  temporal 
variations  in  inlet  fuel  composition  and  temperature. 

The  stack  model  is  a  combination  of  high  order  spatially  dis¬ 
cretized  single  tubular  cell  models.  Compositional  and  temperature 
variations  are  captured  along  each  tube  length.  The  model  captures 
heat  transfer  inside  each  cell,  and  from  a  cell  to  the  gas  outside. 
Thermal  effects  from  cell  to  cell  are  not  included.  To  model  a  stack, 
100  tubular  cells  are  connected  electrically  in  series  and  in  parallel 
with  respect  to  mass  flows.  The  sizes  of  the  tubes  are  taken  as  15  cm 
long  with  an  outer  diameter  of  1  cm.  The  tubular  cell  model  is 
described  in  detail  in  [9], 

The  tail  gas  burner  is  simulated  using  an  axisymmetric  flame 
model  in  Cantera  [10],  The  heat  exchanger  is  modeled  using 
a  dynamic  counter-flow  tubular  heat  exchanger  model  which 
allows  for  varying  temperature  and  composition  in  the  inlet  flows. 
A  detailed  description  of  the  heat  exchanger  model  is  available  in 
[11],  The  physical  parameters  of  the  SOFC  system  are  shown  in 
Table  1.  Together  the  system  model  allows  for  solutions  of  the  gas 
flow  compositions  and  temperatures  throughout  the  system.  A 
diagram  of  the  SOFC  system  and  its  connections  is  shown  in  Fig.  1. 

Analysis  of  the  system  model  provides  excellent  motivation  for 
the  development  of  an  MPC  controller  to  ensure  stable  and  safe 
operation  of  the  system.  For  example,  with  both  the  fuel  reformer 
and  stack  air  flows  supplied  via  blowers,  the  dynamic  response  of 
the  blowers  is  critical  to  the  ability  of  the  reformer  and  stack  to 
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Table  1 

SOFC  system  physical  parameters. 

Parameter  Value 

Fuel  reformer 

Inner  diameter  0.0082  m 

Length  0.0165  m 

SOFC  stack 


Inner  tube  diameter 
Outer  tube  diameter 

Blowers 

Efficiency 


0.15  m 

2205  stainless  steel 
0.25  m 
0.025  m 
0.040  m 
0.005  m 

0.38 

4  x  10  s  kg  m2 


respond  to  demand  changes.  For  large  increases  in  current  demand, 
the  supply  of  fuel  to  the  reformer  can  be  increased  much  more 
quickly  than  the  air,  which  will  result  in  compositions  that 
temporarily  have  very  low  ratios  of  O2/CH4,  and  this  can  lead  to 
carbon  deposition  in  the  stack.  An  increased  current  load  can  also 
produce  very  undesirable  fuel  utilization  because  the  air  blower 
time  constant  is  slower  than  that  of  the  reformer  fuel  supply.  This 
effect  is  shown  in  Fig.  2  for  step  changes  in  between  different 
operating  conditions.  The  simulation  in  Fig.  2  involves  step  changes 
across  a  variety  of  operating  conditions,  from  low  current 
(approximately  5  A)  with  fuel  utilization  varying  from  45%  to  76% 
and  O2/CH4  ratios  from  1.13  to  1.52.  The  “spikes”  present  in  the 
current  and  exhaust  H2%  in  Fig.  2  are  due  to  air  flow  delay  induced 
by  the  blowers.  Large  changes  in  current  require  large  changes  in 
mass  flow  of  both  fuel  and  air.  The  fuel  delivery  system  is  much 
faster  than  the  blowers  delivering  air.  Thus  air/fuel  ratios  can 
become  very  skewed,  resulting  in  undesireable  transients  in  both 
current  and  H2  exhaust.  The  magnitude  of  these  transients  is 
dependent  on  the  dynamics  of  the  blower  model,  but  the  simula¬ 
tion  shows  that  to  avoid  these  transients,  the  controller  must 
account  for  any  mismatch  in  dynamics  between  air  and  fuel 
delivery  subsystems.  An  implemented  controller  must  be  able  to 
choose  between  the  various  mass  flow,  O2/CH4  ratio  combinations 
in  order  to  produce  the  desired  current  without  violating  temper¬ 
ature,  utilization,  and  coking  constraints. 

3.  System  constraints 


control  can  respond  to  current  disturbances  and  maintain  the 
required  stack  temperatures  [12],  However,  compositional 
constraints  have  a  larger  impact  on  the  available  operating  regions. 
These  include  ensuring  the  fuel  utilization,  as  measured  by  the  H2 
percentage  in  the  stack  exhaust,  does  not  decrease  below 
a  minimum  percentage  and  that  the  reformate  composition  is  not 
prone  to  forming  solid  carbon.  These  operating  constraints  are 
summarized  in  Table  2. 

3.1.  Temperature  constraints 

The  CPOX  fuel  reformer  consists  of  a  catalyst  on  a  foam  monolith 
which  encourages  an  exothermic  reaction  in  order  to  produce 
a  suitable  fuel  for  the  stack.  Since  this  reaction  is  exothermic, 
temperatures  can  easily  exceed  the  sintering  temperature  of  the 
catalyst.  Operating  the  reformer  at  too  low  a  temperature,  however, 
can  cause  coking,  which  leads  to  catalyst  deactivation.  Thus  ensuring 
the  temperature  of  the  reformer  is  within  prescribed  limits  forms 
one  of  the  output  constraints  for  the  system  controller.  The 
temperature  limits  for  the  reformer  were  set  to  what  was  experi¬ 
mentally  considered  “safe”  for  a  lab  based  CPOX  reformer  at  the 
Colorado  Fuel  Cell  Center.  The  upper  limit  is  to  avoid  sintering  and 
the  lower  limit  is  the  temperature  the  reformer  was  pre-heated  to 
before  fed  fuel,  in  order  to  alleviate  coking.  In  addition  to  the 
temperature  constraints  on  the  reformer,  the  temperature  of  the 
stack  itself  needs  to  be  regulated  within  a  prescribed  limit.  If  too  low 
a  temperature  is  set,  the  electrolyte  is  no  longer  conductive;  if  the 
stack  temperature  is  too  high,  seals  and  other  components  can  be 
damaged.  The  temperature  limits  for  the  stack  were  set  to  bracket 
the  expected  operating  conditions.  A  lower  limit  is  not  set,  as  the 
controller  ensures  that  the  temperature  is  sufficient  to  produce  the 
demand  current.  The  upper  limit  is  set  to  an  estimated  minimum 
temperature  to  produce  the  operating  range  of  current  (5—25  A). 

3.2.  Hydrogen  utilization 

Fuel  utilization  is  an  important  variable  in  SOFC  system  control. 
Very  low  fuel  utilizations  (and  thus  very  high  exhaust  concentra¬ 
tions  of  H2)  are  to  be  avoided  as  inefficient  and  wasteful.  Very  high 
fuel  utilizations  (very  low  exhaust  H2  concentrations)  indicate  fuel 
depletion  and  can  cause  cell  damage  as  well  as  inability  to  meet 
current  demand.  The  fuel  utilization  is  calculated  as; 

th0utHcomfr  ,out 

titinWcomb  in 


There  are  several  operating  regions  that  need  to  be  avoided  for 
the  SOFC  system.  These  operating  regions  are  dependent  on  both 
the  composition  of  the  gases  and  the  temperatures  of  the  compo¬ 
nents  themselves.  It  is  well  documented  that  model  predictive 


Where  rhout  is  the  mass  flow  of  depleted  fuel,  min  is  the  mass  flow  of 
fuel  into  the  stack  (both  in  kg/sec),  and  Hcomb,x  is  the  heat  of 
combustion  of  whatever  composition  is  being  used.  How  this  utili¬ 
zation  exactly  relates  to  the  H2%  in  the  stack  exhaust  is  a  function  of 
the  composition  of  the  reformate,  the  mass  flows,  and  the  operating 
temperature.  Generally,  however,  if  the  utilization  is  very  high,  then 
the  H2  in  the  stack  will  be  low.  The  primary  control  variable  of 
interest  is  the  H2%  in  the  exhaust.  The  minimum  value  of  2.5%  is 
arbitrarily  chosen  to  prevent  re-oxidation  of  the  anode.  The  actual 
exhaust  H2%  is  allowed  to  vary  in  order  to  meet  transient  current 
demands  as  quickly  as  possible  without  violating  this  minimum.  The 
maximum  H2%  exhaust  concentration  is  not  limited  but  the 
controller  is  designed  to  have  a  preferred  operating  value  of  5%.  This 
exhaust  concentration  corresponds  to  a  fuel  utilization  of  70— 75%. 

3.3.  Reformate  carbon  deposition 

The  deposition  of  solid  carbon  is  analyzed  using  thermodynamic 
equilibrium  chemistry  amongst  the  important  species.  The  method 
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Reformer  Fuel  Mass  Flow  Reformer  Air  Mass  Flow  Stack  Air  Mass  Flow 


Cell  Current  Exhaust  FI  2% 

Fig.  2.  Open  loop  dynamic  response. 


used  is  that  developed  by  [13],  which  was  expanded  to  create 
reforming  fuel  maps  in  [14],  The  determination  of  whether  a  fuel 
can  produce  solid  carbon  is  found  using  the  temperature  and  partial 
pressures  of  H2,  CH4,  CO,  CO2,  and  H2O.  Other  species  are  consid¬ 
ered  inert  and  the  partial  pressures  of  the  active  species  are 
normalized  with  respect  to  the  inert  species.  Considering  the 
following  reactions: 

CH4<-C  +  2H2 

H2  +  2CO<->C  +  H20  (1) 

2CO<->C  +  C02, 

we  calculate  the  conditions  where  the  net  amount  of  C  created  and 
consumed  is  zero.  Graphically,  this  can  be  illustrated  using 
a  ternary  diagram  consisting  of  three  components,  C,  H,  and  0,  as 
shown  in  Fig.  3  (all  thermodynamic  calculations  carried  out  using 
[10]).  The  line  where  the  net  carbon  created  is  zero  is  known  as  the 
carbon  deposition  barrier  (CDB).  It  represents  for  a  given  temper¬ 
ature  the  composition  of  the  reformate  that  is  capable  of  forming 
solid  carbon.  If  the  reformate  composition,  plotted  on  the  same 
ternary  diagram,  is  below  this  line  then  solid  carbon  formation  is 
unlikely.  If  the  reformate  composition  is  above  the  carbon  deposi¬ 
tion  barrier  (CDB)  then  the  reformate  is  prone  to  solid  carbon 
formation.  The  distance  measure  used  is  orthogonal  to  the  H-0 
axis  of  the  ternary  diagram.  This  measure  was  used  as  the  actual 
value  of  interest  is  whether  the  reformate  composition  lies  above  or 
below  the  CDB  curve.  This  analysis  is  based  solely  off  equilibrium 

Table  2 

SOFC  system  constraints. 

Parameter  Value 

MEA  temperature  7/c  <  1200  K 

Reformer  temperature  650  K  <  T„  <  1075  K 

Exhaust  H2%  >2.5% 

CDB  distance  >0.0 


chemistry  and  does  not  consider  the  material  properties  of  the  cell. 
The  provided  measure  of  carbon  formation  is  not  intended  to 
replace  a  detailed  kinematic  analysis.  However,  since  the  object  is 
to  control  the  reformate  composition  it  provides  a  measure  of 
“reformate  quality”  as  far  as  carbon  deposition,  which  has  mostly 
been  ignored  in  previous  SOFC  system  control  methods. 

4.  System  model  reduction 

Model  reduction  via  system  identification  is  a  technique  to  identify 
a  dynamic  model  using  experimental  data.  For  the  previously 
described  SOFC  model,  the  system  is  described  by  a  large  number  of 
non-linear  differential  and  algebraic  differential  equations.  Model 


0  20  40  60  80 

H  O 

Fig.  3.  Fuel  ternary  diagram,  CDB  for  894.76  K  (black),  reformate  composition  with 
inlet  O2/CH4  =  1.97  (red).  (For  interpretation  of  the  references  to  colour  in  this  figure 
legend,  the  reader  is  referred  to  the  web  version  of  this  article.) 
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evaluations  can  take  on  the  order  of  a  minute  for  large  changes  in 
system  state.  Although  this  model  is  easily  accurate  enough  for  control 
purposes,  the  non-linearity  and  evaluation  time  make  it  unsuitable  for 
use  in  a  model  predictive  controller.  Ideally  a  single  linear  model  could 
be  found  that  would  allow  for  a  fast  optimal  control  of  the  system, 
however,  due  to  the  large  non-linearities  present  over  the  operating 
range,  no  such  model  exists.  Thus  the  use  of  a  combination  of  linear 
models  which  are  blended  using  a  function  of  measured  signals  (LPV) 
is  chosen  to  represent  the  SOFC  system. 

Several  approaches  in  the  field  of  system  identification  work  to 
solve  this  problem.  Subspace  LPV  identification  was  chosen  as  it 
creates  a  state  space  model,  suitable  for  model  predictive  control. 
Subspace  LPV  identification  algorithms  can  be  divided  into  two 
approaches,  local  and  global.  Local  approaches  find  linear  models 
around  operating  points  and  interpolate  between  them  using 
a  function  of  a  measured  quantity.  Global  approaches  attempt  to 
excite  all  desired  non-linearities  in  the  system  in  a  single  experi¬ 
ment,  and  then  fit  the  model  based  on  a  functional  dependence  to 
one  or  more  measured  quantities. 

We  use  a  LPV  global  approach  to  the  identification  problem, 
described  in  detail  in  [15]  and  implemented  in  [16].  The  model 
structure  is  given  in  Eq.  (2), 

Xfc+1  =  E  $  (A(i>xk  +  B(‘>uk  +  KUek) 
yk  =  Cxk  +  Duk  +  ek. 

Where  xke  Rn  is  the  state  at  time  k  and  system  matrices  are 
A®e  Rnxn,  BWeRnxr,  JCWeR"**,  Ceilx",  and  DeR,xr  for  a  system 
with  state  size  n,  r  inputs,  and  /  outputs.  /ike  Rmxi  is  a  vector  of 
scheduling  variables  which  determines  how  the  m  different  linear 
models  are  combined  at  each  of  the  L  time  steps.  /if®  is  the  zth 
element  of  /ik.  eke  R;  denotes  the  zero  mean  white  measurement 
residual.  For  the  SOFC  system  analyzed,  model  inputs  are  chosen  as 
the  variables  that  can  be  actuated.  These  include  power  to  the 
reformer  blower  (pbio).  power  to  the  stack  blower  (pstack),  fuel  mass 
flow  to  the  reformer  (mfuei),  and  the  stack  voltage  (V).  Outputs  are 
the  variables  that  require  regulation  or  impose  constraints.  These 
are  taken  as  the  current  (/),  the  membrane  electrode  assembly 
(MEA)  assembly  temperature  (TMea).  the  hydrogen  exhaust 
concentration  (H2,ex)  and  the  temperature  of  the  gas  exiting  the  fuel 
reformer  (Tref)  and  the  distance  from  the  CDB  (dcDB)-  The  input  and 
output  vectors  are  then  defined  as: 

uk  =  [Pbloi  mfuel> Pstack > 

yic  =  [t'mEA-  h  Hi, ex-  Tref.  dcDEi]  • 

The  algorithm  reconstructs  the  state  sequence  using  the  known 
input  and  output  data  as  well  as  the  scheduling  sequence.  Once  this 
is  done  the  unknown  system  matrices  (A®,  B®,  gwt  c,  D,  for 
i  =  {l...m})  are  found  via  least  squares.  The  algorithm  depends  on 
the  existence  of  a  set  of  measurable  scheduling  parameters  (phi), 
which  linearly  weight  the  system  matrices  at  each  time  step.  The 
choice  of  scheduling  parameter  has  a  strong  impact  on  the  accuracy 
of  the  identified  model.  Ideally,  physical  intuition  could  be  used  to 
separate  the  system  into  scheduled  modes  using  a  measured  vari¬ 
able,  however  in  the  case  of  an  SOFC  system  there  are  too  many 
operating  non-linearities  to  make  this  choice  obvious.  The  stack 
current  has  been  used  in  previous  stack  identification  and  control 
applications  [17],  However,  current  does  not  provide  a  perfect 
scheduling  parameter.  There  are  a  multitude  of  operating  points 
with  large  variations  in  reformate  composition,  mass  flow,  and 
voltage,  that  all  produce  the  same  current.  It  is  part  of  the  control 


challenge  to  choose  amongst  these  operating  conditions.  As  such, 
utilizing  only  current  as  a  scheduling  variable  is  limiting,  however, 
additional  measurements  are  often  available.  The  temperature  of 
the  stack  MEA  assembly  can  be  used  as  a  measurement,  as  well  as 
the  temperature  of  the  reformate  gas  leaving  the  fuel  reformer. 
These  additional  measurements  can  provide  additional  compo¬ 
nents  for  the  scheduling  sequence.  Indeed,  previous  work  has 
shown  that  the  inlet  composition  of  the  fuel  to  the  reformer  can  be 
estimated  using  only  the  reformate  temperature  as  a  measurement 
[8[.  The  results  presented  in  this  work  utilize  current  and  MEA 
temperature  as  scheduling  variables. 

For  the  identification  procedure,  variables  to  be  chosen  include 
the  scheduling  parameters,  the  number  of  linear  models  to  blend, 
(m),  the  forward  time  window  (p),  and  the  state  dimension  of  the 
resulting  models,  (n).  Values  selected  are  m  =  3,  p  =  2,  and  n=  7. 
The  model  is  simulated  with  a  time  step  of  1  s.  The  range  of  oper¬ 
ating  conditions  over  which  the  identification  data  is  produced 
must  match  the  desired  operating  regions  of  the  system.  The 
identification  data  set  was  chosen  to  cover  as  wide  a  range  as 
possible:  approximately  0—25  A/cell,  with  both  large  mass  flow  and 
low  mass  flow  conditions  represented,  as  well  as  jumps  between 
them.  Identification  and  validation  simulations  are  plotted  within 
a  reduced  parameter  space  in  Fig.  4.  This  plot  represents  the 
parameter  space  using  three  parameters  which  are  combinations  of 
the  selected  inputs  and  the  output  that  we  desire  to  regulate.  The 
composition  of  fuel  and  air  to  the  reformer  is  shown  using  the 
O2/CH4  ratio.  The  total  mass  flow  of  fuel  and  air  to  both  the  stack 
and  the  blower  and  the  current  are  also  shown.  Using  this 
parameter  space  it  can  be  seen  that  the  identification  simulation  is 
representative  of  the  expected  operating  conditions.  The  identified 
data  points  are  clustered  around  the  steady  state  operating 
conditions  represented  by  the  red  boxes.  If  the  desired  operating 
condition  is  not  close  to  any  of  the  identified  points,  we  would 
expect  the  LPV  model  to  be  a  poor  fit.  This  reduced  parameter  space 
does  not  take  into  account  voltage  or  indicate  which  operating 
conditions  would  violate  system  constraints.  It  does,  however, 
allow  for  a  visualization  of  the  identification  simulation  relative  to 
the  desired  operating  conditions. 

At  each  operating  condition,  pseudo-random  binary  signal 
(PRBS)  perturbations  are  added;  example  inputs  used  are  shown  in 
Fig.  5.  Changes  from  operating  region  to  operating  region  were  also 
simulated,  ensuring  identification  of  the  low  exhaust  H2%  and 
carbon  forming  regions  that  the  controller  needs  to  avoid  during 


Fig.  4.  Identification  and  validation  parameter  space,  the  validation  data  set  is  shown 
using  red  squares.  The  identification  data  set  is  shown  using  blue  circles.  (For  inter¬ 
pretation  of  the  references  to  colour  in  this  figure  legend,  the  reader  is  referred  to  the 
web  version  of  this  article.) 
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actual  operation.  The  scheduling  parameters  (current  and  MEA 
temperature)  are  normalized  to  be  between  negative  and  positive 
one,  and  passed  through  a  low  pass  filter  with  a  normalized  cut  off 
frequency  of  0.4.  The  combination  of  these  measurable  quantities 
allows  for  scheduling  of  an  LPV  model  over  the  operating  range. 
Identification  results  for  a  seven  state,  five  output,  four  input  LPV 
model  are  shown  in  Fig.  6.  The  results  in  Fig.  6  are  given  for  a  vali¬ 
dation  data  set,  using  a  one  step  ahead  predictor.  That  is,  the  state  is 
assumed  to  be  fully  measured  at  each  time  step,  and  the  identified 
model  is  used  to  predict  the  impact  of  the  applied  inputs  to  the 
state  and  the  outputs  at  the  next  time  step.  These  plots  thus  express 
the  accuracy  of  the  LPV  model  when  used  in  a  system  where  all 
outputs  are  measurable  and  the  system  is  fully  observable. 
However,  for  the  SOFC  system  we  are  interested  in  controlling,  only 
three  of  the  outputs  are  measurable.  A  more  useful  representation 
of  the  accuracy  of  the  LPV  model  is  to  combine  the  simulation  with 
a  Kalman  filter,  which  makes  use  of  the  model  and  measurements 
to  estimate  the  unknown  outputs.  The  LPV  model  simulated  with 
a  Kalman  filter  is  shown  in  Fig.  7.  Only  the  two  estimated  outputs 
are  plotted,  as  the  measurements  are  highly  trusted  and  the  model 
is  used  only  to  calculate  the  unknown  outputs. 

The  performance  of  the  estimated  model  is  measured  using 
a  variance  accounted  for  (VAF)  value,  calculated  as: 


f  Mr[yk-yk  ) 

VAFyiI  =  max|l - ^Aoj(lOO).  (4) 

For  the  desired,  4  input,  5  output  system,  identification  of 
a  system  with  three  blended  linear  models  of  state  size  7  resulted  in 
the  highest  VAF  fit  for  a  validation  data  set.  Using  the  current  and 
MEA  temperature  as  scheduling  variables  gives  VAF  scores  (for  the 
one  step  ahead  predictor)  of  97.8927  for  the  MEA  temperature, 


99.5377  for  the  current,  and  99.9619  for  the  reformer  temperature. 
The  remaining  outputs  have  VAF  scores  of  84.2987  for  the  exhaust 
H2%  and  93.8894  for  the  CDB  distance.  These  values  are  for  the 
same  test  run  used  in  the  identification  simulation.  When  testing 
the  identified  model  on  data  sets  other  than  that  used  for  identi¬ 
fication  (but  still  within  the  same  operating  range)  the  VAF 
scores  dropped  for  the  two  unscheduled  variables  to  approximately 
80-90  for  the  exhaust  H2%  and  60  to  90  for  the  CDB  distance.  Both 
temperatures  are  identified  very  well,  with  the  LPV  model 
predicting  the  nonlinear  model  within  0.5  K.  The  current  also 
shows  excellent  identification  with  errors  remaining  within  1  A. 
The  range  of  operation  for  the  identification  simulation  is  also 
quite  challenging,  consisting  of  a  large  number  of  operating 
points  for  the  SOFC  system. 

4.1.  Model  predictive  control 

Model  predictive  control  has  been  applied  to  SOFC  systems 
previously  [12,18,19]  but  not  with  a  model  of  sufficient  fidelity  to 
estimate  fuel  composition  throughout.  The  use  of  an  LPV  model 
allows  for  linear  MPC  to  be  applied  where  previous  work  has 
resorted  to  nonlinear  MPC  implementations.  With  the  development 
of  a  sufficiently  detailed  system  level  model,  carbon  formation  can 
be  avoided  and  temperature  limits  maintained  while  still  allowing 
load  following.  The  controller  signal  flow  diagram  is  shown  in  Fig.  8. 
The  measurable  outputs  are  taken  from  the  SOFC  system  at  each 
time  step.  The  scheduling  sequence  is  calculated  from  the  measured 
outputs  and  used  to  calculate  the  current  linear  model  using  the  LPV 
framework.  This  linear  model  is  then  used  within  a  standard  Kal¬ 
man  filter  (along  with  the  measurements)  to  obtain  a  full  state 
estimate.  This  gives  an  estimate  for  the  two  unmeasurable  quanti¬ 
ties  we  are  interested  in  (H2  exhaust  and  the  CDB  distance).  With  an 
estimate  for  these  two  quantities  in  hand,  the  model  at  the  current 
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Fig.  5.  Inputs  used  for  identification  simulation. 
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Fig.  6.  System  identification  results,  identified  model  is  solid  and  black,  nonlinear  model  is  red  and  dashed.  (For  interpretation  of  the  references  to  colour  in  this  figure  legend,  the 
reader  is  referred  to  the  web  version  of  this  article.) 


time-step  calculated  and  the  desired  current,  the  standard  MPC 
optimization  problem  can  be  solved  over  a  time  window  of  N 
samples.  The  resulting  quadratic  cost  optimization  problem  (Eq.  (5)) 
gives  the  inputs  to  be  applied  at  the  next  time-step. 


Note  that  variables  are  as  defined  in  Eq  (2).  A  tracking  imple¬ 
mentation  with  a  N  =  7  step  horizon  is  implemented  with  output 
weights  only  on  the  current,  exhaust  H2%  and  the  CDB  distance. 
The  demand  current,  target  exhaust  H2%  and  target  CDB  distance 
are  considered  static  over  the  horizon  (1...N)  for  the  optimization 


J  =  minimize  x'kNPfxkN  +  £  xTklQxk  l  +  uTkiRukyi  +  \\Qy (. Vkj  ~ yref,i)  1 1 2 

subject  to  yki(l)  <  1200, i  =  0,  ...,N  -  1 

ykJ(  3)  >  .025,i  bO . N-l 

650  <yfci(4)  <  1075, i  =  0,...,N-  1 
yM(5)>0,i  =  0,...,JV-l 


a 


b 


Time  (s) 


Time  (s) 
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Fig.  7.  LPV  model  results  with  Kalman  filter.  The  LPV  and  Kalman  estimate  is  black  ar 
colour  in  this  figure  legend,  the  reader  is  referred  to  the  web  version  of  this  article.) 
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Fig.  8.  MPC  process  flow. 


problem  in  Eq.  (5).  They  are  also  updated  at  every  time  step  after 
applying  the  calculated  inputs.  The  weighting  for  the  current  is 
five  orders  of  magnitude  larger  than  that  of  the  exhaust  H2%  and 
CDB  distance.  This  weighting  has  the  effect  of  primarily  ensuring 
the  demand  current  is  met.  The  smaller  weights  push  the 
controller  to  keep  the  exhaust  H2%  and  CDB  distance  above  their 
hard  constraint  minimums.  This  stabilization  is  to  account  for 
modeling  error,  forcing  the  controller  to  operate  at  some  distance 
from  the  hard  constraints.  The  state  and  input  weights  are  both  set 
to  identity  as  the  only  requirement  is  to  stay  within  hard 
constraints.  In  order  to  enforce  stability,  a  terminal  weight  is 
added  to  the  states.  This  weighting  matrix  Pf  is  chosen  as  a  solu¬ 
tion  to  a  Riccati  inequality  [20],  All  weights  are  summarized  in 
Table  3. 


Value 


Qy 


R 


Q 


0  0  0  0  0 

0  to5  0  0  0 

0  0  1.0  0  0 

0  0  0  0  0 

0  0  0  0  1.0 

'1  0  0  0] 

0  10  0 

0  0  10 
0  0  0  lj 
"1  0  0  0  0 
0  10  0  0 

0  0  10  0 
0  0  0  1  0 
0  0  0  0  1 


5.  Results 

The  controller  was  tested  over  a  variety  of  output  conditions. 
Results  are  shown  for  a  current  demand  step  both  down  and  up  of 
1.0  A  per  second  and  0.1  A  per  second  in  Fig.  9.  In  both  cases  the 
controller  is  able  to  closely  meet  the  current  demand,  however  the 
fast  current  step  causes  a  much  larger  oscillation  (Fig.  9(a))  than  for 
the  slow  demand  change  (Fig.  9(d)).  The  output  constraints  are  also 
violated  for  the  large  demand  change,  shown  in  Fig.  9(b)  and  (c). 
Violations  occur  for  the  exhaust  H 2%  constraint  during  the  large 
current  demand  decrease.  While  violations  of  the  CDB  constraint 
occur  for  both  current  demand  increases  and  decreases.  However, 
the  controller  can  stabilize  the  output  current  over  the  entire 
operating  range  (5-25  A).  The  controller  is  thus  capable  of  load 
following  while  not  violating  any  of  the  system  constraints,  for 
a  limited  rate  of  current  change.  The  reason  that  the  current 
trajectory  rate  change  must  be  limited  is  that  the  MPC  controller 
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Fig.  10.  k  =  7  step  predictor,  identified  model  is  solid  and  black,  nonlinear  model  is  red  and  dashed.  (For  interpretation  of  the  references  to  colour  in  this  figure  legend,  the  reader  is 
referred  to  the  web  version  of  this  article.) 


operates  over  a  prediction  horizon.  The  model  used  for  the 
optimization  is  generated  from  the  LPV  model  and  fixed  over  this 
prediction  horizon.  A  fixed  linear  model  greatly  decreases  the 
complexity  of  the  resulting  optimization  problem,  but  introduces 
modeling  error.  To  visualize  the  model  discrepencies  faced  by  the 
controller,  we  utilize  a  k-step  predictor  where  k  is  equal  to  the 
MPC  horizon.  The  results  of  a  k  =  7  step  predictor  are  shown  in 
Fig.  10.  The  controller  model  faces  large  transient  errors  in  H2 
exhaust  concentration  and  CDB  distance  over  the  horizon,  as 
shown  in  Fig.  10(c)  and  (e).  These  results  are  for  the  same  input 
sequence  as  that  used  in  Figs.  6  and  7. 

The  MPC  controller  thus  has  a  discrepency  between  the  pre¬ 
dicted  dynamics,  which  are  used  for  the  optimization,  and  the 
dynamics  which  are  used  for  the  optimization  at  each  successive 
time  step.  If  the  current  demanded  is  allowed  to  undergo  large  step 
changes,  the  dynamics  at  the  next  time-step  may  be  different  than 
predicted,  and  thus  the  calculated  inputs  will  show  a  large  error. 
The  controller  will  eventually  converge  to  a  new  model,  but  oscil¬ 
lation  will  occur  around  the  new  set-points.  The  magnitude  and 
length  of  the  oscillation  is  dependent  on  how  far  away  the  current 
MPC  model  is  and  thus  how  inaccurate  the  effects  of  the  calculated 
outputs.  Violation  of  hard  constraints  can  occur  during  these 
changes  as  well  for  the  same  reasons.  The  MPC  still  calculates 
inputs  that  do  not  cause  output  violations,  but  the  model  it  uses  to 
do  so  is  inaccurate  over  the  horizon.  In  addition  this  mismatch  is 
not  equal  at  all  operating  points.  For  some  regions,  the  LPV  model  is 
very  close,  thus  allowing  fast  current  changes,  in  other  regions,  the 
LPV  model  is  inaccurate,  requiring  a  limitation  on  the  allowed 
magnitude  of  current  change.  The  effect  of  this  mismatch  is 
prominently  displayed  in  errors  between  the  MPC  predicted 
current  (or  the  current  demanded)  and  the  actual  system  current 
which  is  measured. 

Minimizing  the  MPC  prediction  horizon  reduces  the  effect  of 
this  error,  but  increases  the  controller  costs  and  can  result  in  an 
infeasible  problem.  Implementing  some  form  of  LPV  MPC  should 


greatly  reduce  this  problem,  although  with  a  cost  of  increasing  the 
complexity  of  the  controller.  Alternatively,  the  rate  of  change  for 
the  scheduling  variables  can  be  limited.  For  the  SOFC  system,  this 
requires  limiting  the  allowed  rate  of  current  demand  change.  The 
controller  has  no  difficulty  avoiding  constraints  with  a  current 
demand  change  of  0.1  A.  The  use  of  a  battery  or  capacitor  will  be 
required  to  limit  the  magnitude  of  the  desired  load  change  to  such 
a  value. 

6.  Conclusion 

Current  challenges  in  SOFC  operation  are  often  related  to 
durability  and  carbon  formation,  this  work  has  demonstrated 
a  controller  which  can  mitigate  these  design  challenges.  The 
controller  is  capable  of  staying  within  all  desired  operating 
temperatures  as  well  as  ensuring  the  reformer  produces  a  fuel 
which  is  not  prone  to  carbon  formation,  provided  the  magnitude  of 
the  load  change  is  limited.  The  system  identification  algorithm 
requires  only  an  existing  first  principles  model  which  can  produce 
input/output  data  over  the  operating  range.  Linear  MPC  is  used, 
which  is  both  fast  and  well  established. 
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